#include "MagicFormula.h"
#include <math.h>

MagicFormula::MagicFormula(){
}

MagicFormula::~MagicFormula(){
}

void MagicFormula::InitialPrivateFx(double parameters[]){
    parameters[0] = 1.5;
    parameters[1] = 1;
    parameters[2] = 0;
    parameters[3] = 1;
    parameters[4] = -1;
    parameters[5] = 0.5;
    parameters[6] = 0;
    parameters[7] = 0.1;
    parameters[8] = 15;
    parameters[9] = 0;
    parameters[10] = 0;
    parameters[11] = 0;
    parameters[12] = 0;
    parameters[13] = 0;
    parameters[14] = 0;
}

void MagicFormula::InitialPrivateFy(double parameters[]){
    parameters[0] =1.3507;
    parameters[1] =1;
    parameters[2] =-0.18;
    parameters[3] =-2.88;
    parameters[4] =-0.0074;
    parameters[5] =-0.0063;
    parameters[6] =-0.5;
    parameters[7] =0;
    parameters[8] =-15;
    parameters[9] =2;
    parameters[10] =-0.025;
    parameters[11] =0;
    parameters[12] =0;
    parameters[13] =0;
    parameters[14] =0;
    parameters[15] =0;
    parameters[16] =0;
    parameters[17] =0;
    /*
    parameters[0] = +1.362068411445;
    parameters[1] = +0.861928153584;
    parameters[2] = -0.218584421478;
    parameters[3] = +0.291165288703;
    parameters[4] = -4.474358292612;
    parameters[5] = +4.462274617102;
    parameters[6] = -0.077786812099;
    parameters[7] = +3.745114441482;
    parameters[8] = -6.335948555537;
    parameters[9] = +1.020567635342;
    parameters[10] = +0.190148868942;
    parameters[11] = -0.003046960023;
    parameters[12] = -0.006145878982;
    parameters[13] = -0.074795516942;
    parameters[14] = +0.002608298293;
    parameters[15] = -0.007846849532;
    parameters[16] = +0.357479162362;
    parameters[17] = +0.042236285149;
    */
}

void MagicFormula::InitialPrivateMz(double parameters[]){
    //★其他的系数都是用ComputePrivateFy()求好的，在这个函数里都是常数。所以这个函数共有25个系数
    parameters[0] = 10.904;
    parameters[1] = -1.8412;
    parameters[2] = -0.52041;
    parameters[3] = 0.03921;
    parameters[4] = 0.41511;
    parameters[5] = 8.9846;
    parameters[6] = 0;
    parameters[7] = 1.2136;
    parameters[8] = 1.093509;
    parameters[9] = -0.0092183;
    parameters[10] = -0.057061;
    parameters[11] = 0.73954;
    parameters[12] = -0.0067783;
    parameters[13] = 0.0052254;
    parameters[14] = -0.18175;
    parameters[15] = 0.029952;
    parameters[16] = -1.5697;
    parameters[17] = -0.33394;
    parameters[18] = 0;
    parameters[19] = 0.26711;
    parameters[20] = -3.594;
    parameters[21] = 0.0047326;
    parameters[22] = 0.0026687;
    parameters[23] = 0.11998;
    parameters[24] = 0.059083;
    /*
    parameters[0]= 8.934839027233;
    parameters[1]= -6.941364385261;
    parameters[2]= 2.104262608423;
    parameters[3]= 0.126403102967;
    parameters[4]= 0.001816792888;
    parameters[5]= 10.545361460877;
    parameters[6]= -0.043727090934;
    parameters[7]= 1.113596518029;
    parameters[8]= 0.167646503391;
    parameters[9]= 0.074444156645;
    parameters[10]= -0.130776460886;
    parameters[11]= -1.243282315173;
    parameters[12]= 0.011015950328;
    parameters[13]= 0.005993392241;
    parameters[14]= 0.192996932579;
    parameters[15]= 0.156353001746;
    parameters[16]= -2.557048457299;
    parameters[17]= -1.319917490676;
    parameters[18]= -0.985950424383;
    parameters[19]= -0.719406526948;
    parameters[20]= -4.580848256598;
    parameters[21]= 0.012793420027;
    parameters[22]= -0.000661038371;
    parameters[23]= 0.041855143008;
    parameters[24]= -0.029155456862;
    */
}

void MagicFormula::InitialPrivateCombFx(double parameters[]){
    parameters[0] = 13.276;
    parameters[1] = -13.778;
    parameters[2] = 1.2568;
    parameters[3] = 0.65225;
    parameters[4] = -0.24948;
    parameters[5] = 0.0050722;
}

void MagicFormula::InitialPrivateCombFy(double parameters[]){
    parameters[0] = 7.1433;
    parameters[1] = 9.1916;
    parameters[2] = -0.027856;
    parameters[3] = 1.0719;
    parameters[4] = -0.27572;
    parameters[5] = 0.32802;
    parameters[6] = 5.7448e-006;
    parameters[7] = -3.1368e-005;
    parameters[8] = -0.027825;
    parameters[9] = 0.053604;
    parameters[10] = -0.27568;
    parameters[11] = 12.12;
    parameters[12] = 1.9;
    parameters[13] = -10.704;
}

void MagicFormula::InitialPrivateCombMz(double parameters[]){
    parameters[1] = 0.033372;
    parameters[2] = 0.0043624;
    parameters[3] = 0.56742;
    parameters[4] = -0.24116;
}
//以下是便于从mtalab移植方程而设置的宏
#define Fzo dataPrivate[0]
#define dfz dataPrivate[1]
//以下是PrivateFx的中间变量
#define Shx dataPrivateFx[0]
#define Svx dataPrivateFx[1]
#define kx dataPrivateFx[2]
#define Cx dataPrivateFx[3]
#define Dx dataPrivateFx[4]
#define Ex dataPrivateFx[5]
#define CFx dataPrivateFx[6]
#define Bx dataPrivateFx[7]
#define Fxo dataPrivateFx[8]
//以下是PrivateFy的中间变量
#define Shy dataPrivateFy[0]
#define Svy dataPrivateFy[1]
#define alpha_y dataPrivateFy[2]
#define Cy dataPrivateFy[3]
#define uy dataPrivateFy[4]
#define Dy dataPrivateFy[5]
#define Ey dataPrivateFy[6]
#define CFy dataPrivateFy[7]
#define By dataPrivateFy[8]
#define Fyo dataPrivateFy[9]
//PrivateMz的中间变量，PI_2 = 2/pi
#define PI_2 0.63661977236758134307553505349006
#define Ro dataPrivateMz[0]
#define alpha_t dataPrivateMz[1]
#define Bt dataPrivateMz[2]
#define Ct dataPrivateMz[3]
#define Dt dataPrivateMz[4]
#define Et dataPrivateMz[5]
#define t dataPrivateMz[6]
#define alpha_r dataPrivateMz[7]
#define Br dataPrivateMz[8]
#define Cr dataPrivateMz[9]
#define Dr dataPrivateMz[10]
#define Mzr dataPrivateMz[11]
//PrivateCombFx函数的中间变量
#define alpha_s dataPrivateCombFx[0]
#define Bxa dataPrivateCombFx[1]
#define Exa dataPrivateCombFx[2]
#define Gxa1 dataPrivateCombFx[3]
#define Gxa2 dataPrivateCombFx[4]
#define CombFxsim0 dataPrivateCombFx[5]
//PrivateCombFy函数的中间变量
#define Shyk dataPrivateCombFy[0]
#define k_s dataPrivateCombFy[1]
#define Byk dataPrivateCombFy[2]
#define Eyk dataPrivateCombFy[3]
#define Svyk dataPrivateCombFy[4]
#define Gyk1 dataPrivateCombFy[5]
#define Gyk2 dataPrivateCombFy[6]
#define CombFysim0 dataPrivateCombFy[7]
//PrivateCombMz函数的中间变量
#define alpha_teq dataPrivateCombMz[0]
#define tMz dataPrivateCombMz[1]
#define alpha_req dataPrivateCombMz[2]
#define MzrMz dataPrivateCombMz[3]
#define s dataPrivateCombMz[4]

//PrivateFx函数和其偏导函数
double MagicFormula::PrivateFx(double Fz, double k, double r, double parameters[]){
    /*与matlab源代码的系数对应关系表:
    #define Pcx1 parameters[0]
    #define Pdx1 parameters[1]
    #define Pdx2 parameters[2]
    #define Pdx3 parameters[3]
    #define Pex1 parameters[4]
    #define Pex2 parameters[5]
    #define Pex3 parameters[6]
    #define Pex4 parameters[7]
    #define Pkx1 parameters[8]
    #define Pkx2 parameters[9]
    #define Pkx3 parameters[10]
    #define Phx1 parameters[11]
    #define Phx2 parameters[12]
    #define Pvx1 parameters[13]
    #define Pvx2 parameters[14]
    */
    /*matlab源程序方程代码:
    Fzo=11500;
    dfz=(Fz-Fzo)./Fzo;
    Shx=(Phx1+Phx2*dfz);
    Svx=Fz.*(Pvx1+Pvx2*dfz);
    kx=k+Shx;
    Cx=Pcx1;
    ux=(Pdx1+Pdx2*dfz).*(1-Pdx3*r.^2);
    Dx=ux.*Fz;
    Ex=(Pex1+Pex2*dfz+Pex3*dfz.^2).*(1-Pex4*sign(kx));
    CFk=Fz.*(Pkx1+Pkx2*dfz).*exp(Pkx3*dfz);
    Bx=CFk./(Cx.*Dx);
    Fxo  = (Dx.*sin(Cx.*atan(Bx.*kx-Ex.*(Bx.*kx-atan(Bx.*kx)))))+Svx;
    Fxsim=Fxo;
    */
    //    Fzo=11500;	//Fzo提前输入进去
    dfz = Fz / Fzo - 1;
    Shx = (parameters[11] + parameters[12] * dfz);
    Svx = Fz*(parameters[13] + parameters[14] * dfz);
    kx = k + Shx;
    Dx = Fz*(parameters[1] + parameters[2] * dfz)*(1 - parameters[3] * r*r);
    Ex = (parameters[4] + parameters[5] * dfz + parameters[6] * dfz*dfz)*(1 - parameters[7] * sign(kx));
    CFx = Fz*(parameters[8] + parameters[9] * dfz)*exp(parameters[10] * dfz);
    Bx = CFx / (parameters[0] * Dx);
    Fxo = Dx*sin(parameters[0] * atan(Bx*kx - Ex*(Bx*kx - atan(Bx*kx)))) + Svx;
    return Fxo;
}

double MagicFormula::PrivateFx_00(double Fz, double k, double r, double parameters[]){//checked by matlab
    const double Bx_ = -CFx*Dx / ((parameters[0] * Dx)*(parameters[0] * Dx));
    const double temp1 = Bx*kx - Ex*(Bx*kx - atan(Bx*kx));
    const double temp1_ = kx*Bx_ - Ex*kx*Bx_ + Ex / (1 + Bx*Bx*kx*kx)*kx*Bx_;
    const double temp2 = atan(temp1);
    const double temp2_ = 1 / (1 + temp1*temp1)*temp1_;
    return Dx*cos(parameters[0] * temp2)*(temp2 + parameters[0] * temp2_);	//=Fxo_  for parameters[0]
}

double MagicFormula::PrivateFx_01(double Fz, double k, double r, double parameters[]){//checked by matlab
    const double Dx_ = Fz*(1 - parameters[3] * r*r);
    const double Bx_ = -CFx*parameters[0] / ((parameters[0] * Dx)*(parameters[0] * Dx))*Dx_;
    const double temp1 = Bx*kx - Ex*(Bx*kx - atan(Bx*kx));
    const double temp1_ = kx*Bx_ - Ex*kx*Bx_ + Ex / (1 + Bx*Bx*kx*kx) *kx*Bx_;
    const double temp2 = sin(parameters[0] * atan(temp1));
    const double temp2_ = cos(parameters[0] * atan(temp1)) *parameters[0] / (1 + temp1*temp1) *temp1_;
    return Dx_*temp2 + temp2_*Dx;		//=Fxo_  for parameters[1]
}

double MagicFormula::PrivateFx_02(double Fz, double k, double r, double parameters[]){//checked by matlab
    const double Dx_ = Fz*dfz*(1 - parameters[3] * r*r);
    const double Bx_ = -CFx*parameters[0] / ((parameters[0] * Dx)*(parameters[0] * Dx))*Dx_;
    const double temp1 = Bx*kx - Ex*(Bx*kx - atan(Bx*kx));
    const double temp1_ = kx*Bx_ - Ex*kx*Bx_ + Ex / (1 + Bx*Bx*kx*kx) *kx*Bx_;
    const double temp2 = sin(parameters[0] * atan(temp1));
    const double temp2_ = cos(parameters[0] * atan(temp1)) *parameters[0] / (1 + temp1*temp1) *temp1_;
    return Dx_*temp2 + temp2_*Dx;		//=Fxo_  for parameters[2]
}

double MagicFormula::PrivateFx_03(double Fz, double k, double r, double parameters[]){//checked by matlab
    const double Dx_ = -Fz*(parameters[1] + parameters[2] * dfz)*r*r;
    const double Bx_ = -dfz*parameters[0] / ((parameters[0] * Dx)*(parameters[0] * Dx))*Dx_;
    const double temp1 = Bx*kx - Ex*(Bx*kx - atan(Bx*kx));
    const double temp1_ = kx*Bx_ - Ex*kx*Bx_ + Ex / (1 + Bx*Bx*kx*kx) *kx*Bx_;
    const double temp2 = sin(parameters[0] * atan(temp1));
    const double temp2_ = cos(parameters[0] * atan(temp1)) *parameters[0] / (1 + temp1*temp1) *temp1_;
    return Dx_*temp2 + temp2_*Dx;		//=Fxo_  for parameters[3]
}

double MagicFormula::PrivateFx_04(double Fz, double k, double r, double parameters[]){//checked by matlab
    const double Ex_ = (1 - parameters[7] * sign(kx));
    const double temp1 = Bx*kx - Ex*(Bx*kx - atan(Bx*kx));
    const double temp1_ = -Ex_*(Bx*kx - atan(Bx*kx));
    return Dx*cos(parameters[0] * atan(temp1)) *parameters[0] / (1 + temp1*temp1) *temp1_;//=Fxo_  for parameters[4]
}

double MagicFormula::PrivateFx_05(double Fz, double k, double r, double parameters[]){//checked by matlab
    const double Ex_ = dfz*(1 - parameters[7] * sign(kx));
    const double temp1 = Bx*kx - Ex*(Bx*kx - atan(Bx*kx));
    const double temp1_ = -Ex_*(Bx*kx - atan(Bx*kx));
    return Dx*cos(parameters[0] * atan(temp1)) *parameters[0] / (1 + temp1*temp1) *temp1_;//=Fxo_  for parameters[5]
}

double MagicFormula::PrivateFx_06(double Fz, double k, double r, double parameters[]){//checked by matlab
    const double Ex_ = dfz*dfz*(1 - parameters[7] * sign(kx));
    const double temp1 = Bx*kx - Ex*(Bx*kx - atan(Bx*kx));
    const double temp1_ = -Ex_*(Bx*kx - atan(Bx*kx));
    return Dx*cos(parameters[0] * atan(temp1)) *parameters[0] / (1 + temp1*temp1) *temp1_;//=Fxo_  for parameters[6]
}

double MagicFormula::PrivateFx_07(double Fz, double k, double r, double parameters[]){//checked by matlab
    const double Ex_ = -sign(kx)*(parameters[4] + parameters[5] * dfz + parameters[6] * dfz*dfz);
    const double temp1 = Bx*kx - Ex*(Bx*kx - atan(Bx*kx));
    const double temp1_ = -Ex_*(Bx*kx - atan(Bx*kx));
    return Dx*cos(parameters[0] * atan(temp1)) *parameters[0] / (1 + temp1*temp1) *temp1_;//=Fxo_  for parameters[7]
}

double MagicFormula::PrivateFx_08(double Fz, double k, double r, double parameters[]){//checked by matlab
    const double CFx_ = Fz*exp(parameters[10] * dfz);
    const double Bx_ = CFx_ / (parameters[0] * Dx);
    const double temp1 = Bx*kx - Ex*(Bx*kx - atan(Bx*kx));
    const double temp1_ = kx*Bx_ - Ex*kx*Bx_ + Ex / (1 + Bx*Bx*kx*kx) *kx*Bx_;
    return Dx*cos(parameters[0] * atan(temp1)) *parameters[0] / (1 + temp1*temp1) *temp1_;//=Fxo_  for parameters[8]
}

double MagicFormula::PrivateFx_09(double Fz, double k, double r, double parameters[]){//checked by matlab
    const double CFx_ = Fz*dfz*exp(parameters[10] * dfz);
    const double Bx_ = CFx_ / (parameters[0] * Dx);
    const double temp1 = Bx*kx - Ex*(Bx*kx - atan(Bx*kx));
    const double temp1_ = kx*Bx_ - Ex*kx*Bx_ + Ex / (1 + Bx*Bx*kx*kx) *kx*Bx_;
    return Dx*cos(parameters[0] * atan(temp1)) *parameters[0] / (1 + temp1*temp1) *temp1_;//=Fxo_  for parameters[9]
}

double MagicFormula::PrivateFx_10(double Fz, double k, double r, double parameters[]){//checked by matlab
    const double CFx_ = Fz*(parameters[8] + parameters[9] * dfz)*exp(parameters[10] * dfz)*dfz;
    const double Bx_ = CFx_ / (parameters[0] * Dx);
    const double temp1 = Bx*kx - Ex*(Bx*kx - atan(Bx*kx));
    const double temp1_ = kx*Bx_ - Ex*kx*Bx_ + Ex / (1 + Bx*kx*Bx*kx) *kx*Bx_;
    return Dx*cos(parameters[0] * atan(temp1)) *parameters[0] / (1 + temp1*temp1) *temp1_;//=Fxo_  for parameters[10]
}

double MagicFormula::PrivateFx_11(double Fz, double k, double r, double parameters[]){//checked by matlab
    const double kx_ = 1;
    const double temp1 = Bx*kx - Ex*(Bx*kx - atan(Bx*kx));
    const double temp1_ = kx_*Bx - Ex*kx_*Bx + Ex / (1 + Bx*Bx*kx*kx) *kx_*Bx;
    return Dx*cos(parameters[0] * atan(temp1)) *parameters[0] / (1 + temp1*temp1) *temp1_;//=Fxo_  for parameters[11]
}

double MagicFormula::PrivateFx_12(double Fz, double k, double r, double parameters[]){//checked by matlab
    const double kx_ = dfz;
    const double temp1 = Bx*kx - Ex*(Bx*kx - atan(Bx*kx));
    const double temp1_ = kx_*Bx - Ex*kx_*Bx + Ex / (1 + Bx*kx*Bx*kx) *kx_*Bx;
    return Dx*cos(parameters[0] * atan(temp1)) *parameters[0] / (1 + temp1*temp1) *temp1_;//=Fxo_  for parameters[12]
}

double MagicFormula::PrivateFx_13(double Fz, double k, double r, double parameters[]){//checked by matlab
    return Fz;		//=Fxo_  for parameters[13]
}

double MagicFormula::PrivateFx_14(double Fz, double k, double r, double parameters[]){	//checked by matlab
    return Fz*dfz;		//=Fxo_  for parameters[14]
}

double MagicFormula::PrivateFy(double Fz, double alpha, double r, double parameters[]){//checked by matlab
    /*与matlab源代码的系数对应关系表:
    Pcy1=parametersFy[0];
    Pdy1=parametersFy[1];
    Pdy2=parametersFy[2];
    Pdy3=parametersFy[3];
    Pey1=parametersFy[4];
    Pey2=parametersFy[5];
    Pey3=parametersFy[6];
    Pey4=parametersFy[7];
    Pky1=parametersFy[8];
    Pky2=parametersFy[9];
    Pky3=parametersFy[10];
    Phy1=parametersFy[11];
    Phy2=parametersFy[12];
    Phy3=parametersFy[13];
    Pvy1=parametersFy[14];
    Pvy2=parametersFy[15];
    Pvy3=parametersFy[16];
    Pvy4=parametersFy[17];
    */
    /*matlab源程序方程代码:
    Fzo=11500;
    dfz=(Fz-Fzo)./Fzo;
    Shy=(Phy1+Phy2*dfz)+Phy3*r;
    Svy=Fz.*((Pvy1+Pvy2*dfz)+(Pvy3+Pvy4*dfz).*r);
    alpha_y=alpha+Shy;
    Cy=Pcy1;
    uy=(Pdy1+Pdy2*dfz).*(1-Pdy3*r.^2);
    Dy=uy.*Fz;
    Ey=(Pey1+Pey2*dfz).*(1-(Pey3+Pey4*r).*sign(alpha_y));
    CFyo=Pky1*Fzo.*sin(2*atan(Fz./(Pky2*Fzo)));
    CFy=CFyo.*(1-Pky3*abs(r));
    By=CFy./(Cy.*Dy);
    Fyo  = (Dy.*sin(Cy.*atan(By.*alpha_y-Ey.*(By.*alpha_y-atan(By.*alpha_y)))))+Svy;
    Fysim=Fyo;
    */
    //	Fzo=11500;	//Fzo提前输入进去
    dfz = Fz / Fzo - 1;
    Shy = (parameters[11] + parameters[12] * dfz) + parameters[13] * r;
    Svy = Fz*((parameters[14] + parameters[15] * dfz) + (parameters[16] + parameters[17] * dfz)*r);
    alpha_y = alpha + Shy;
    Cy = parameters[0];
    uy = (parameters[1] + parameters[2] * dfz)*(1 - parameters[3] * r*r);
    Dy = Fz*uy;
    Ey = (parameters[4] + parameters[5] * dfz)*(1 - (parameters[6] + parameters[7] * r)*sign(alpha_y));
    CFy = parameters[8] * Fzo*sin(2 * atan(Fz / (parameters[9] * Fzo)))*(1 - parameters[10] * fabs(r));
    By = CFy / (parameters[0] * Dy);
    Fyo = (Dy*sin(parameters[0] * atan(By*alpha_y - Ey*(By*alpha_y - atan(By*alpha_y))))) + Svy;
    return Fyo;
}

double MagicFormula::PrivateFy_00(double parameters[]){//checked by matlab
    const double By_ = -CFy*Dy / (parameters[0] * Dy*parameters[0] * Dy);
    const double temp1 = By*alpha_y - Ey*(By*alpha_y - atan(By*alpha_y));
    const double temp1_ = alpha_y*By_ - alpha_y*Ey*By_ + Ey / (1 + By*alpha_y*By*alpha_y)*alpha_y*By_;
    const double temp2 = atan(temp1);
    const double temp2_ = 1 / (1 + temp1*temp1)*temp1_;
    return Dy*cos(parameters[0] * temp2)*(temp2 + parameters[0] * temp2_);
}

double MagicFormula::PrivateFy_01(double Fz, double r, double parameters[]){//checked by matlab
    const double Dy_ = Fz*(1 - parameters[3] * r*r);
    const double By_ = -CFy*parameters[0] / ((parameters[0] * Dy)*(parameters[0] * Dy))*Dy_;
    const double temp1 = By*alpha_y - Ey*(By*alpha_y - atan(By*alpha_y));
    const double temp1_ = alpha_y*By_ - alpha_y*Ey*By_ + Ey / (1 + By*alpha_y*By*alpha_y)*alpha_y*By_;
    const double temp2 = sin(parameters[0] * atan(temp1));
    const double temp2_ = cos(parameters[0] * atan(temp1)) *parameters[0] / (1 + temp1*temp1) *temp1_;
    return Dy_*temp2 + temp2_*Dy;
}

double MagicFormula::PrivateFy_02(double Fz, double r, double parameters[]){//checked by matlab
    const double Dy_ = Fz*dfz*(1 - parameters[3] * r*r);
    const double By_ = -CFy*parameters[0] / ((parameters[0] * Dy)*(parameters[0] * Dy))*Dy_;
    const double temp1 = By*alpha_y - Ey*(By*alpha_y - atan(By*alpha_y));
    const double temp1_ = alpha_y*By_ - alpha_y*Ey*By_ + Ey / (1 + By*alpha_y*By*alpha_y)*alpha_y*By_;
    const double temp2 = sin(parameters[0] * atan(temp1));
    const double temp2_ = cos(parameters[0] * atan(temp1)) *parameters[0] / (1 + temp1*temp1) *temp1_;
    return Dy_*temp2 + temp2_*Dy;
}

double MagicFormula::PrivateFy_03(double Fz, double r, double parameters[]){//checked by matlab
    const double Dy_ = -Fz*(parameters[1] + parameters[2] * dfz)*r*r;
    const double By_ = -CFy*parameters[0] / ((parameters[0] * Dy)*(parameters[0] * Dy))*Dy_;
    const double temp1 = By*alpha_y - Ey*(By*alpha_y - atan(By*alpha_y));
    const double temp1_ = alpha_y*By_ - alpha_y*Ey*By_ + Ey / (1 + By*alpha_y*By*alpha_y)*alpha_y*By_;
    const double temp2 = sin(parameters[0] * atan(temp1));
    const double temp2_ = cos(parameters[0] * atan(temp1)) *parameters[0] / (1 + temp1*temp1) *temp1_;
    return Dy_*temp2 + temp2_*Dy;
}

double MagicFormula::PrivateFy_04(double r, double parameters[]){//checked by matlab
    const double Ey_ = 1 - (parameters[6] + parameters[7] * r)*sign(alpha_y);
    const double temp1 = By*alpha_y - Ey*(By*alpha_y - atan(By*alpha_y));
    const double temp1_ = -Ey_*(By*alpha_y - atan(By*alpha_y));
    return Dy*cos(parameters[0] * atan(temp1)) *parameters[0] / (1 + temp1*temp1) *temp1_;
}

double MagicFormula::PrivateFy_05(double r, double parameters[]){//checked by matlab
    const double Ey_ = dfz*(1 - (parameters[6] + parameters[7] * r)*sign(alpha_y));
    const double temp1 = By*alpha_y - Ey*(By*alpha_y - atan(By*alpha_y));
    const double temp1_ = -Ey_*(By*alpha_y - atan(By*alpha_y));
    return Dy*cos(parameters[0] * atan(temp1)) *parameters[0] / (1 + temp1*temp1) *temp1_;
}

double MagicFormula::PrivateFy_06(double parameters[]){//error！corrected
    const double Ey_ = -(parameters[4] + parameters[5] * dfz)*sign(alpha_y);
    const double temp1 = By*alpha_y - Ey*(By*alpha_y - atan(By*alpha_y));
    const double temp1_ = -Ey_*(By*alpha_y - atan(By*alpha_y));
    return Dy*cos(parameters[0] * atan(temp1)) *parameters[0] / (1 + temp1*temp1) *temp1_;
}

double MagicFormula::PrivateFy_07(double r, double parameters[]){//error！corrected
    const double Ey_ = -(parameters[4] + parameters[5] * dfz)*r*sign(alpha_y);
    const double temp1 = By*alpha_y - Ey*(By*alpha_y - atan(By*alpha_y));
    const double temp1_ = -Ey_*(By*alpha_y - atan(By*alpha_y));
    return Dy*cos(parameters[0] * atan(temp1)) *parameters[0] / (1 + temp1*temp1) *temp1_;
}

double MagicFormula::PrivateFy_08(double Fz, double r, double parameters[]){//checked by matlab
    const double CFy_ = Fzo*sin(2 * atan(Fz / (parameters[9] * Fzo)))*(1 - parameters[10] * fabs(r));
    const double By_ = CFy_ / (parameters[0] * Dy);
    const double temp1 = By*alpha_y - Ey*(By*alpha_y - atan(By*alpha_y));
    const double temp1_ = alpha_y*By_ - alpha_y*Ey*By_ + Ey / (1 + By*alpha_y*By*alpha_y)*alpha_y*By_;
    return Dy*cos(parameters[0] * atan(temp1)) *parameters[0] / (1 + temp1*temp1) *temp1_;
}

double MagicFormula::PrivateFy_09(double Fz, double r, double parameters[]){//error! corrected
    const double temp0 = Fz / (parameters[9] * Fzo);
    const double temp0_ = -Fz / (parameters[9] * Fzo*parameters[9] * Fzo)*Fzo;
    const double CFy_ = parameters[8] * Fzo*(1 - parameters[10] * fabs(r))*cos(2 * atan(temp0)) * 2 / (1 + temp0*temp0)*temp0_;
    const double By_ = CFy_ / (parameters[0] * Dy);
    const double temp1 = By*alpha_y - Ey*(By*alpha_y - atan(By*alpha_y));
    const double temp1_ = alpha_y*By_ - alpha_y*Ey*By_ + Ey / (1 + By*alpha_y*By*alpha_y)*alpha_y*By_;
    return Dy*cos(parameters[0] * atan(temp1)) *parameters[0] / (1 + temp1*temp1) *temp1_;
}

double MagicFormula::PrivateFy_10(double Fz, double r, double parameters[]){//checked by matlab
    const double CFy_ = -parameters[8] * Fzo*sin(2 * atan(Fz / (parameters[9] * Fzo)))*fabs(r);
    const double By_ = CFy_ / (parameters[0] * Dy);
    const double temp1 = By*alpha_y - Ey*(By*alpha_y - atan(By*alpha_y));
    const double temp1_ = alpha_y*By_ - alpha_y*Ey*By_ + Ey / (1 + By*alpha_y*By*alpha_y)*alpha_y*By_;
    return Dy*cos(parameters[0] * atan(temp1)) *parameters[0] / (1 + temp1*temp1) *temp1_;
}

double MagicFormula::PrivateFy_11(double parameters[]){//checked by matlab
    const double alpha_y_ = 1;
    const double temp1 = By*alpha_y - Ey*(By*alpha_y - atan(By*alpha_y));
    const double temp1_ = By*alpha_y_ - Ey*By*alpha_y_ + Ey / (1 + By*alpha_y*By*alpha_y)*By*alpha_y_;
    return Dy*cos(parameters[0] * atan(temp1)) *parameters[0] / (1 + temp1*temp1) *temp1_;
}

double MagicFormula::PrivateFy_12(double parameters[]){//checked by matlab
    const double alpha_y_ = dfz;
    const double temp1 = By*alpha_y - Ey*(By*alpha_y - atan(By*alpha_y));
    const double temp1_ = By*alpha_y_ - Ey*By*alpha_y_ + Ey / (1 + By*alpha_y*By*alpha_y)*By*alpha_y_;
    return Dy*cos(parameters[0] * atan(temp1)) *parameters[0] / (1 + temp1*temp1) *temp1_;
}

double MagicFormula::PrivateFy_13(double r, double parameters[]){//checked by matlab
    const double alpha_y_ = r;
    const double temp1 = By*alpha_y - Ey*(By*alpha_y - atan(By*alpha_y));
    const double temp1_ = By*alpha_y_ - Ey*By*alpha_y_ + Ey / (1 + By*alpha_y*By*alpha_y)*By*alpha_y_;
    return Dy*cos(parameters[0] * atan(temp1)) *parameters[0] / (1 + temp1*temp1) *temp1_;
}

double MagicFormula::PrivateFy_14(double Fz){
    return Fz;
}

double MagicFormula::PrivateFy_15(double Fz){
    return Fz*dfz;
}

double MagicFormula::PrivateFy_16(double Fz, double r){
    return Fz*r;
}

double MagicFormula::PrivateFy_17(double Fz, double r){
    return Fz*r*dfz;
}


double MagicFormula::PrivateMz(double Fz, double alpha, double r, double parameters[]){//checked with matlab
    /*与matlab源代码的系数对应关系表:
    qbz1=parameters[0];
    qbz2=parameters[1];
    qbz3=parameters[2];
    qbz4=parameters[3];
    qbz5=parameters[4];
    qbz9=parameters[5];
    qbz10=parameters[6];
    qcz1=parameters[7];
    qdz1=parameters[8];
    qdz2=parameters[9];
    qdz3=parameters[10];
    qdz4=parameters[11];
    qdz6=parameters[12];
    qdz7=parameters[13];
    qdz8=parameters[14];
    qdz9=parameters[15];
    qez1=parameters[16];
    qez2=parameters[17];
    qez3=parameters[18];
    qez4=parameters[19];
    qez5=parameters[20];
    qhz1=parameters[21];
    qhz2=parameters[22];
    qhz3=parameters[23];
    qhz4=parameters[24];
    */
    /*matlab源程序方程代码:
    Fzo=11500;
    dfz=(Fz-Fzo)./Fzo;
    Shy=(Phy1+Phy2*dfz)+Phy3*r;
    Svy=Fz.*((Pvy1+Pvy2*dfz)+(Pvy3+Pvy4*dfz).*r);
    alpha_y=alpha+Shy;
    Cy=Pcy1;
    uy=(Pdy1+Pdy2*dfz).*(1-Pdy3*r.^2);
    Dy=uy.*Fz;
    Ey=(Pey1+Pey2*dfz).*(1-(Pey3+Pey4*r).*sign(alpha_y));
    CFyo=Pky1*Fzo.*sin(2*atan(Fz./(Pky2*Fzo)));
    CFy=CFyo.*(1-Pky3*abs(r));
    By=CFy./(Cy.*Dy);
    Fyo  = (Dy.*sin(Cy.*atan(By.*alpha_y-Ey.*(By.*alpha_y-atan(By.*alpha_y)))))+Svy;
    %上面和PrivateFy一样
    Sht=qhz1+qhz2*dfz+(qhz3+qhz4*dfz).*r;
    alpha_t=tan(alpha)+Sht;
    Bt=(qbz1+qbz2*dfz+qbz3*dfz.^2).*(1+qbz4*r+qbz5*abs(r));
    Ct=qcz1;
    Dt=Fz.*(qdz1+qdz2*dfz).*(1+qdz3*r+qdz4*r.^2).*(Ro./Fzo);
    Et=(qez1+qez2*dfz+qez3*dfz.^2).*(1+(qez4+qez5*r).*((2/pi)*atan(Bt.*Ct.*alpha_t)));
    t=Dt.*cos(Ct.*atan(Bt.*alpha_t-Et.*(Bt.*alpha_t-atan(Bt.*alpha_t)))).*cos(alpha);
    Shf=Shy+Svy./CFy;
    alpha_r=tan(alpha)+Shf;
    Br=qbz9+qbz10.*By.*Cy;
    Cr=1;
    Dr=Fz.*((qdz6+qdz7*dfz)+(qdz8+qdz9*dfz).*r).*Ro;
    Mzr=Dr.*cos(Cr.*atan(Br.*alpha_r)).*cos(alpha);
    Mzo=-t.*Fyo+Mzr;
    Mzsim=Mzo;
    */
    /*
    const double Fzo=11500;
    const double dfz=(Fz-Fzo)/Fzo;
    const double Shy=(parametersFy[11]+parametersFy[12]*dfz)+parametersFy[13]*r;
    const double Svy=Fz*((parametersFy[14]+parametersFy[15]*dfz)+(parametersFy[16]+parametersFy[17]*dfz)*r);
    const double alpha_y=alpha+Shy;
    const double Dy=Fz*(parametersFy[1]+parametersFy[2]*dfz)*(1-parametersFy[3]*r*r);
    const double Ey=(parametersFy[4]+parametersFy[5]*dfz)*(1-(parametersFy[6]+parametersFy[7]*r)*sign(alpha_y));
    const double CFy=parametersFy[8]*Fzo*sin(2*atan(Fz/(parametersFy[9]*Fzo)))*(1-parametersFy[10]*fabs(r));
    const double By=CFy/(parametersFy[0]*Dy);
    const double Fyo=(Dy*sin(parametersFy[0]*atan(By*alpha_y-Ey*(By*alpha_y-atan(By*alpha_y)))))+Svy;

    const double alpha_t=tan(alpha)+qhz1+qhz2*dfz+(qhz3+qhz4*dfz)*r;
    const double Bt=(qbz1+qbz2*dfz+qbz3*dfz*dfz)*(1+qbz4*r+qbz5*fabs(r));
    const double Dt=Fz*(qdz1+qdz2*dfz)*(1+qdz3*r+qdz4*r*r)*(Ro/Fzo);
    const double Et=(qez1+qez2*dfz+qez3*dfz*dfz)*(1+(qez4+qez5*r)*((2/pi)*atan(Bt*qcz1*alpha_t)));
    const double t=Dt*cos(qcz1*atan(Bt*alpha_t-Et*(Bt*alpha_t-atan(Bt*alpha_t))))*cos(alpha);
    const double alpha_r=tan(alpha)+Shy+Svy/CFy;
    const double Br=qbz9+qbz10*By*parametersFy[0];
    const double Cr=1;
    const double Ro=0.375;
    const double Dr=Fz*((qdz6+qdz7*dfz)+(qdz8+qdz9*dfz)*r)*Ro;
    const double Mzr=Dr*cos(Cr*atan(Br*alpha_r))*cos(alpha);
    return =-t*Fyo+Mzr;
    */

    PrivateFy(Fz, alpha, r, parametersPrivateFy);//计算PrivateFy用到的中间数据
    Cr = 1;
    alpha_t = tan(alpha) + parameters[21] + parameters[22] * dfz + (parameters[23] + parameters[24] * dfz)*r;
    Bt = (parameters[0] + parameters[1] * dfz + parameters[2] * dfz*dfz)*(1 + parameters[3] * r + parameters[4] * fabs(r));
    Ct = parameters[7];
    Dt = Fz*(parameters[8] + parameters[9] * dfz)*(1 + parameters[10] * r + parameters[11] * r*r)*(Ro / Fzo);
    Et = (parameters[16] + parameters[17] * dfz + parameters[18] * dfz*dfz)*(1 + (parameters[19] + parameters[20] * r)*(PI_2*atan(Bt*parameters[7] * alpha_t)));
    t = Dt*cos(parameters[7] * atan(Bt*alpha_t - Et*(Bt*alpha_t - atan(Bt*alpha_t))))*cos(alpha);
    alpha_r = tan(alpha) + Shy + Svy / CFy;
    Br = parameters[5] + parameters[6] * By*Cy;
    Dr = Fz*Ro*((parameters[12] + parameters[13] * dfz) + (parameters[14] + parameters[15] * dfz)*r);
    Mzr = Dr*cos(Cr*atan(Br*alpha_r))*cos(alpha);
    return -t*Fyo + Mzr;
}

double MagicFormula::PrivateMz_00(double alpha, double r, double parameters[]){//只精确到4位有效数字
    const double Bt_ = (1 + parameters[3] * r + parameters[4] * fabs(r));
    const double temp1 = Bt*alpha_t - Et*(Bt*alpha_t - atan(Bt*alpha_t));
    const double temp1_ = alpha_t*Bt_ - Et*alpha_t*Bt_ + Et / (1 + Bt*alpha_t*Bt*alpha_t)*alpha_t*Bt_;
    const double t_ = -Dt*cos(alpha)*sin(parameters[7] * atan(temp1))*parameters[7] / (1 + temp1*temp1)*temp1_;
    return -t_*Fyo;
}

double MagicFormula::PrivateMz_01(double alpha, double r, double parameters[]){//只精确到4位有效数字
    const double Bt_ = dfz*(1 + parameters[3] * r + parameters[4] * fabs(r));
    const double temp1 = Bt*alpha_t - Et*(Bt*alpha_t - atan(Bt*alpha_t));
    const double temp1_ = alpha_t*Bt_ - Et*alpha_t*Bt_ + Et / (1 + Bt*alpha_t*Bt*alpha_t)*alpha_t*Bt_;
    const double t_ = -Dt*cos(alpha)*sin(parameters[7] * atan(temp1))*parameters[7] / (1 + temp1*temp1)*temp1_;
    return -t_*Fyo;
}

double MagicFormula::PrivateMz_02(double alpha, double r, double parameters[]){//只精确到4位有效数字
    const double Bt_ = dfz*dfz*(1 + parameters[3] * r + parameters[4] * fabs(r));
    const double temp1 = Bt*alpha_t - Et*(Bt*alpha_t - atan(Bt*alpha_t));
    const double temp1_ = alpha_t*Bt_ - Et*alpha_t*Bt_ + Et / (1 + Bt*alpha_t*Bt*alpha_t)*alpha_t*Bt_;
    const double t_ = -Dt*cos(alpha)*sin(parameters[7] * atan(temp1))*parameters[7] / (1 + temp1*temp1)*temp1_;
    return -t_*Fyo;
}

double MagicFormula::PrivateMz_03(double alpha, double r, double parameters[]){//只精确到4位有效数字
    const double Bt_ = r*(parameters[0] + parameters[1] * dfz + parameters[2] * dfz*dfz);
    const double temp1 = Bt*alpha_t - Et*(Bt*alpha_t - atan(Bt*alpha_t));
    const double temp1_ = alpha_t*Bt_ - Et*alpha_t*Bt_ + Et / (1 + Bt*alpha_t*Bt*alpha_t)*alpha_t*Bt_;
    const double t_ = -Dt*cos(alpha)*sin(parameters[7] * atan(temp1))*parameters[7] / (1 + temp1*temp1)*temp1_;
    return -t_*Fyo;
}

double MagicFormula::PrivateMz_04(double alpha, double r, double parameters[]){//只精确到4位有效数字
    const double Bt_ = fabs(r)*(parameters[0] + parameters[1] * dfz + parameters[2] * dfz*dfz);
    const double temp1 = Bt*alpha_t - Et*(Bt*alpha_t - atan(Bt*alpha_t));
    const double temp1_ = alpha_t*Bt_ - Et*alpha_t*Bt_ + Et / (1 + Bt*alpha_t*Bt*alpha_t)*alpha_t*Bt_;
    const double t_ = -Dt*cos(alpha)*sin(parameters[7] * atan(temp1))*parameters[7] / (1 + temp1*temp1)*temp1_;
    return -t_*Fyo;
}

double MagicFormula::PrivateMz_05(double alpha){//checked by matlab
    const double Br_ = 1;
    return -Dr*sin(Cr*atan(Br*alpha_r))*cos(alpha)*Cr / (1 + Br*alpha_r*Br*alpha_r)*alpha_r*Br_;
}

double MagicFormula::PrivateMz_06(double alpha){//checked by matlab
    const double Br_ = By*Cy;
    return -Dr*sin(Cr*atan(Br*alpha_r))*cos(alpha)*Cr / (1 + Br*alpha_r*Br*alpha_r)*alpha_r*Br_;
}

double MagicFormula::PrivateMz_07(double alpha, double r, double parameters[]){//只精确到9位有效数字
    const double Et_ = (parameters[16] + parameters[17] * dfz + parameters[18] * dfz*dfz)*(parameters[19] + parameters[20] * r)*PI_2*alpha_t / (1 + Bt*parameters[7] * Bt*parameters[7])*Bt;
    const double temp1 = Bt*alpha_t - Et*(Bt*alpha_t - atan(Bt*alpha_t));
    const double temp1_ = -(Bt*alpha_t - atan(Bt*alpha_t))*Et_;
    const double t_ = -Dt*sin(parameters[7] * atan(temp1))*cos(alpha)*(atan(temp1) + parameters[7] / (1 + temp1*temp1)*temp1_);
    return -t_*Fyo;
}

double MagicFormula::PrivateMz_08(double Fz, double alpha, double r, double parameters[]){//checked by matlab
    const double Dt_ = Fz*(1 + parameters[10] * r + parameters[11] * r*r)*(Ro / Fzo);
    const double t_ = Dt_*cos(parameters[7] * atan(Bt*alpha_t - Et*(Bt*alpha_t - atan(Bt*alpha_t))))*cos(alpha);
    return -t_*Fyo;
}

double MagicFormula::PrivateMz_09(double Fz, double alpha, double r, double parameters[]){//checked by matlab
    const double Dt_ = Fz*dfz*(1 + parameters[10] * r + parameters[11] * r*r)*(Ro / Fzo);
    const double t_ = Dt_*cos(parameters[7] * atan(Bt*alpha_t - Et*(Bt*alpha_t - atan(Bt*alpha_t))))*cos(alpha);
    return -t_*Fyo;
}

double MagicFormula::PrivateMz_10(double Fz, double alpha, double r, double parameters[]){//checked by matlab
    const double Dt_ = Fz*(parameters[8] + parameters[9] * dfz)*r*(Ro / Fzo);
    const double t_ = Dt_*cos(parameters[7] * atan(Bt*alpha_t - Et*(Bt*alpha_t - atan(Bt*alpha_t))))*cos(alpha);
    return -t_*Fyo;
}

double MagicFormula::PrivateMz_11(double Fz, double alpha, double r, double parameters[]){//checked by matlab
    const double Dt_ = Fz*(parameters[8] + parameters[9] * dfz)*r*r*(Ro / Fzo);
    const double t_ = Dt_*cos(parameters[7] * atan(Bt*alpha_t - Et*(Bt*alpha_t - atan(Bt*alpha_t))))*cos(alpha);
    return -t_*Fyo;
}

double MagicFormula::PrivateMz_12(double Fz, double alpha){
    const double Dr_ = Fz*Ro;
    return Dr_*cos(Cr*atan(Br*alpha_r))*cos(alpha);
}

double MagicFormula::PrivateMz_13(double Fz, double alpha){//checked by matlab
    const double Dr_ = Fz*Ro*dfz;
    return Dr_*cos(Cr*atan(Br*alpha_r))*cos(alpha);
}

double MagicFormula::PrivateMz_14(double Fz, double alpha, double r){//checked by matlab
    const double Dr_ = Fz*Ro*r;
    return Dr_*cos(Cr*atan(Br*alpha_r))*cos(alpha);
}

double MagicFormula::PrivateMz_15(double Fz, double alpha, double r){//checked by matlab
    const double Dr_ = Fz*Ro*dfz*r;
    return Dr_*cos(Cr*atan(Br*alpha_r))*cos(alpha);
}

double MagicFormula::PrivateMz_16(double alpha, double r, double parameters[]){//error! corrected
    const double Et_ = 1 + (parameters[19] + parameters[20] * r)*(PI_2*atan(Bt*parameters[7] * alpha_t));
    const double temp1 = Bt*alpha_t - Et*(Bt*alpha_t - atan(Bt*alpha_t));
    const double temp1_ = -Et_*(Bt*alpha_t - atan(Bt*alpha_t));
    const double t_ = -Dt*cos(alpha)*sin(parameters[7] * atan(temp1))*parameters[7] / (1 + temp1*temp1)*temp1_;
    return -t_*Fyo;
}

double MagicFormula::PrivateMz_17(double alpha, double r, double parameters[]){//error! corrected
    const double Et_ = dfz*(1 + (parameters[19] + parameters[20] * r)*(PI_2*atan(Bt*parameters[7] * alpha_t)));
    const double temp1 = Bt*alpha_t - Et*(Bt*alpha_t - atan(Bt*alpha_t));
    const double temp1_ = -Et_*(Bt*alpha_t - atan(Bt*alpha_t));
    const double t_ = -Dt*sin(parameters[7] * atan(temp1))*cos(alpha)*parameters[7] / (1 + temp1*temp1)*temp1_;
    return -t_*Fyo;
}

double MagicFormula::PrivateMz_18(double alpha, double r, double parameters[]){//error! corrected
    const double Et_ = dfz*dfz*(1 + (parameters[19] + parameters[20] * r)*(PI_2*atan(Bt*parameters[7] * alpha_t)));
    const double temp1 = Bt*alpha_t - Et*(Bt*alpha_t - atan(Bt*alpha_t));
    const double temp1_ = -Et_*(Bt*alpha_t - atan(Bt*alpha_t));
    const double t_ = -Dt*sin(parameters[7] * atan(temp1))*cos(alpha)*parameters[7] / (1 + temp1*temp1)*temp1_;
    return -t_*Fyo;
}

double MagicFormula::PrivateMz_19(double alpha, double parameters[]){//error! corrected
    const double Et_ = (parameters[16] + parameters[17] * dfz + parameters[18] * dfz*dfz)*(PI_2*atan(Bt*parameters[7] * alpha_t));
    const double temp1 = Bt*alpha_t - Et*(Bt*alpha_t - atan(Bt*alpha_t));
    const double temp1_ = -Et_*(Bt*alpha_t - atan(Bt*alpha_t));
    const double t_ = -Dt*sin(parameters[7] * atan(temp1))*cos(alpha)*parameters[7] / (1 + temp1*temp1)*temp1_;
    return -t_*Fyo;
}

double MagicFormula::PrivateMz_20(double alpha, double r, double parameters[]){//error! corrected
    const double Et_ = (parameters[16] + parameters[17] * dfz + parameters[18] * dfz*dfz)*r*(PI_2*atan(Bt*parameters[7] * alpha_t));
    const double temp1 = Bt*alpha_t - Et*(Bt*alpha_t - atan(Bt*alpha_t));
    const double temp1_ = -Et_*(Bt*alpha_t - atan(Bt*alpha_t));
    const double t_ = -Dt*sin(parameters[7] * atan(temp1))*cos(alpha)*parameters[7] / (1 + temp1*temp1)*temp1_;
    return -t_*Fyo;
}

double MagicFormula::PrivateMz_21(double alpha, double parameters[]){//只精确到4位有效数字
    const double alpha_t_ = 1;
    const double temp1 = Bt*alpha_t - Et*(Bt*alpha_t - atan(Bt*alpha_t));
    const double temp1_ = Bt*alpha_t_ - Et*Bt*alpha_t_ + Et / (1 + Bt*alpha_t*Bt*alpha_t)*Bt*alpha_t_;
    const double t_ = -Dt*sin(parameters[7] * atan(temp1))*cos(alpha)*parameters[7] / (1 + temp1*temp1)*temp1_;
    return -t_*Fyo;
}

double MagicFormula::PrivateMz_22(double alpha, double parameters[]){//只精确到4位有效数字
    const double alpha_t_ = dfz;
    const double temp1 = Bt*alpha_t - Et*(Bt*alpha_t - atan(Bt*alpha_t));
    const double temp1_ = Bt*alpha_t_ - Et*Bt*alpha_t_ + Et / (1 + Bt*alpha_t*Bt*alpha_t)*Bt*alpha_t_;
    const double t_ = -Dt*sin(parameters[7] * atan(temp1))*cos(alpha)*parameters[7] / (1 + temp1*temp1)*temp1_;
    return -t_*Fyo;
}

double MagicFormula::PrivateMz_23(double alpha, double r, double parameters[]){//只精确到4位有效数字
    const double alpha_t_ = r;
    const double temp1 = Bt*alpha_t - Et*(Bt*alpha_t - atan(Bt*alpha_t));
    const double temp1_ = Bt*alpha_t_ - Et*Bt*alpha_t_ + Et / (1 + Bt*alpha_t*Bt*alpha_t)*Bt*alpha_t_;
    const double t_ = -Dt*sin(parameters[7] * atan(temp1))*cos(alpha)*parameters[7] / (1 + temp1*temp1)*temp1_;
    return -t_*Fyo;
}

double MagicFormula::PrivateMz_24(double alpha, double r, double parameters[]){//只精确到4位有效数字
    const double alpha_t_ = dfz*r;
    const double temp1 = Bt*alpha_t - Et*(Bt*alpha_t - atan(Bt*alpha_t));
    const double temp1_ = Bt*alpha_t_ - Et*Bt*alpha_t_ + Et / (1 + Bt*alpha_t*Bt*alpha_t)*Bt*alpha_t_;
    const double t_ = -Dt*sin(parameters[7] * atan(temp1))*cos(alpha)*parameters[7] / (1 + temp1*temp1)*temp1_;
    return -t_*Fyo;
}

double MagicFormula::PrivateCombFx(double Fz, double k, double alpha, double r, double parameters[]){
    /*与matlab源代码的系数对应关系表:
    Rbx1=parameters[0];
    Rbx2=parameters[1];
    Rcx1=parameters[2];
    Rex1=parameters[3];
    Rex2=parameters[4];
    Rhx1=parameters[5];
    */
    /*matlab源程序方程代码:
    Fzo=4000;
    dfz=(Fz-Fzo)./Fzo;
    Shxa=Rhx1;
    alpha_s=alpha+Shxa;
    Bxa=Rbx1*cos(atan(Rbx2*k));
    Cxa=Rcx1;
    Exa=Rex1+Rex2*dfz;
    Gxa=cos(Cxa*atan(Bxa.*alpha_s-Exa.*(Bxa.*alpha_s-atan(Bxa.*alpha_s))))./cos(Cxa*atan(Bxa.*Shxa-Exa.*(Bxa.*Shxa-atan(Bxa.*Shxa))));
    Fxo=PrivateFx(Fz,k,r,parametersFx);
    Fxsim=Fxo.*Gxa;
    */
    //先用最优参数把PrivateFx的Fxo算出来
    PrivateFx(Fz, k, r, parametersPrivateFx);	//计算Fxo中间变量
    alpha_s = alpha + parameters[5];
    Bxa = parameters[0] * cos(atan(parameters[1] * k));
    Exa = parameters[3] + parameters[4] * dfz;
    Gxa1 = cos(parameters[2] * atan(Bxa*alpha_s - Exa*(Bxa*alpha_s - atan(Bxa*alpha_s))));
    Gxa2 = cos(parameters[2] * atan(Bxa*parameters[5] - Exa*(Bxa*parameters[5] - atan(Bxa*parameters[5]))));
    CombFxsim0 = Fxo*Gxa1 / Gxa2;
    return CombFxsim0;
}

double MagicFormula::PrivateCombFx_00(double k, double parameters[]){//checked by matlab
    const double Bxa_ = cos(atan(parameters[1] * k));
    const double temp1 = Bxa*alpha_s - Exa*(Bxa*alpha_s - atan(Bxa*alpha_s));
    const double temp1_ = alpha_s*Bxa_ - Exa*alpha_s*Bxa_ + Exa / (1 + Bxa*alpha_s*Bxa*alpha_s)*alpha_s*Bxa_;
    const double Gxa1_ = -sin(parameters[2] * atan(temp1))*parameters[2] / (1 + temp1*temp1)*temp1_;
    const double temp2 = Bxa*parameters[5] - Exa*(Bxa*parameters[5] - atan(Bxa*parameters[5]));
    const double temp2_ = parameters[5] * Bxa_ - Exa*parameters[5] * Bxa_ + Exa / (1 + Bxa*parameters[5] * Bxa*parameters[5])*parameters[5] * Bxa_;
    const double Gxa2_ = -sin(parameters[2] * atan(temp2))*parameters[2] / (1 + temp2*temp2)*temp2_;
    return Fxo*(Gxa1_*Gxa2 - Gxa1*Gxa2_) / (Gxa2*Gxa2);
}

double MagicFormula::PrivateCombFx_01(double k, double parameters[]){//checked by matlab
    const double Bxa_ = -parameters[0] * sin(atan(parameters[1] * k)) / (1 + parameters[1] * k*parameters[1] * k)*k;
    const double temp1 = Bxa*alpha_s - Exa*(Bxa*alpha_s - atan(Bxa*alpha_s));
    const double temp1_ = alpha_s*Bxa_ - Exa*alpha_s*Bxa_ + Exa / (1 + Bxa*alpha_s*Bxa*alpha_s)*alpha_s*Bxa_;
    const double Gxa1_ = -sin(parameters[2] * atan(temp1))*parameters[2] / (1 + temp1*temp1)*temp1_;
    const double temp2 = Bxa*parameters[5] - Exa*(Bxa*parameters[5] - atan(Bxa*parameters[5]));
    const double temp2_ = parameters[5] * Bxa_ - Exa*parameters[5] * Bxa_ + Exa / (1 + Bxa*parameters[5] * Bxa*parameters[5])*parameters[5] * Bxa_;
    const double Gxa2_ = -sin(parameters[2] * atan(temp2))*parameters[2] / (1 + temp2*temp2)*temp2_;
    return Fxo*(Gxa1_*Gxa2 - Gxa1*Gxa2_) / (Gxa2*Gxa2);
}

double MagicFormula::PrivateCombFx_02(double parameters[]){//checked by matlab
    const double temp1 = atan(Bxa*alpha_s - Exa*(Bxa*alpha_s - atan(Bxa*alpha_s)));
    const double Gxa1_ = -sin(parameters[2] * temp1)*temp1;
    const double temp2 = atan(Bxa*parameters[5] - Exa*(Bxa*parameters[5] - atan(Bxa*parameters[5])));
    const double Gxa2_ = -sin(parameters[2] * temp2)*temp2;
    return Fxo*(Gxa1_*Gxa2 - Gxa1*Gxa2_) / (Gxa2*Gxa2);
}

double MagicFormula::PrivateCombFx_03(double parameters[]){//checked by matlab
    const double Exa_ = 1;
    const double temp1 = Bxa*alpha_s - Exa*(Bxa*alpha_s - atan(Bxa*alpha_s));
    const double temp1_ = -(Bxa*alpha_s - atan(Bxa*alpha_s))*Exa_;
    const double Gxa1_ = -sin(parameters[2] * atan(temp1))*parameters[2] / (1 + temp1*temp1)*temp1_;
    const double temp2 = Bxa*parameters[5] - Exa*(Bxa*parameters[5] - atan(Bxa*parameters[5]));
    const double temp2_ = -(Bxa*parameters[5] - atan(Bxa*parameters[5]))*Exa_;
    const double Gxa2_ = -sin(parameters[2] * atan(temp2))*parameters[2] / (1 + temp2*temp2)*temp2_;
    return Fxo*(Gxa1_*Gxa2 - Gxa1*Gxa2_) / (Gxa2*Gxa2);
}

double MagicFormula::PrivateCombFx_04(double parameters[]){//checked by matlab
    const double Exa_ = dfz;
    const double temp1 = Bxa*alpha_s - Exa*(Bxa*alpha_s - atan(Bxa*alpha_s));
    const double temp1_ = -(Bxa*alpha_s - atan(Bxa*alpha_s))*Exa_;
    const double Gxa1_ = -sin(parameters[2] * atan(temp1))*parameters[2] / (1 + temp1*temp1)*temp1_;
    const double temp2 = Bxa*parameters[5] - Exa*(Bxa*parameters[5] - atan(Bxa*parameters[5]));
    const double temp2_ = -(Bxa*parameters[5] - atan(Bxa*parameters[5]))*Exa_;
    const double Gxa2_ = -sin(parameters[2] * atan(temp2))*parameters[2] / (1 + temp2*temp2)*temp2_;
    return Fxo*(Gxa1_*Gxa2 - Gxa1*Gxa2_) / (Gxa2*Gxa2);
}

double MagicFormula::PrivateCombFx_05(double parameters[]){//checked by matlab
    const double alpha_s_ = 1;
    const double temp1 = Bxa*alpha_s - Exa*(Bxa*alpha_s - atan(Bxa*alpha_s));
    const double temp1_ = Bxa*alpha_s_ - Exa*Bxa*alpha_s_ + Exa / (1 + Bxa*alpha_s*Bxa*alpha_s)*Bxa*alpha_s_;
    const double Gxa1_ = -sin(parameters[2] * atan(temp1))*parameters[2] / (1 + temp1*temp1)*temp1_;
    const double temp2 = Bxa*parameters[5] - Exa*(Bxa*parameters[5] - atan(Bxa*parameters[5]));
    const double temp2_ = Bxa - Exa*Bxa + Exa / (1 + Bxa*parameters[5] * Bxa*parameters[5])*Bxa;
    const double Gxa2_ = -sin(parameters[2] * atan(temp2))*parameters[2] / (1 + temp2*temp2)*temp2_;
    return Fxo*(Gxa1_*Gxa2 - Gxa1*Gxa2_) / (Gxa2*Gxa2);
}

double MagicFormula::PrivateCombFy(double Fz, double k, double alpha, double r, double parameters[]){
    /*与matlab源代码的系数对应关系表:
    Rby1=parameters[0];
    Rby2=parameters[1];
    Rby3=parameters[2];
    Rcy1=parameters[3];
    Rey1=parameters[4];
    Rey2=parameters[5];
    Rhy1=parameters[6];
    Rhy2=parameters[7];
    Rvy1=parameters[8];
    Rvy2=parameters[9];
    Rvy3=parameters[10];
    Rvy4=parameters[11];
    Rvy5=parameters[12];
    Rvy6=parameters[13];
    */
    /*matlab源程序方程代码:
    Fzo=9179;
    dfz=(Fz-Fzo)./Fzo;
    uy=(Pdy1+Pdy2*dfz).*(1-Pdy3*r.^2);
    Shyk=Rhy1+Rhy2.*dfz;
    k_s=k+Shyk;
    Byk=Rby1.*cos(atan(Rby2.*(alpha-Rby3)));
    Cyk=Rcy1;
    Eyk=Rey1+Rey2.*dfz;
    Dvyk=uy.*Fz.*(Rvy1+Rvy2.*dfz+Rvy3.*r).*cos(atan(Rvy4.*alpha));
    Svyk=Dvyk.*sin(Rvy5.*atan(Rvy6.*k));
    Gyk=cos(Cyk.*atan(Byk.*k_s-Eyk.*(Byk.*k_s-atan(Byk.*k_s))))./cos(Cyk.*atan(Byk.*Shyk-Eyk.*(Byk.*Shyk-atan(Byk.*Shyk))));
    Fyo=PrivateFy(Fz,alpha,r,parametersFy);
    Fysim=Fyo.*Gyk+Svyk;
    */
    PrivateFy(Fz, alpha, r, parametersPrivateFy);	//计算uy、Fyo两个中间变量
    Shyk = parameters[6] + parameters[7] * dfz;
    k_s = k + Shyk;
    Byk = parameters[0] * cos(atan(parameters[1] * (alpha - parameters[2])));
    Eyk = parameters[4] + parameters[5] * dfz;
    Svyk = uy*Fz*(parameters[8] + parameters[9] * dfz + parameters[10] * r)*cos(atan(parameters[11] * alpha))*sin(parameters[12] * atan(parameters[13] * k));
    Gyk1 = cos(parameters[3] * atan(Byk*k_s - Eyk*(Byk*k_s - atan(Byk*k_s))));
    Gyk2 = cos(parameters[3] * atan(Byk*Shyk - Eyk*(Byk*Shyk - atan(Byk*Shyk))));
    CombFysim0 = Fyo*Gyk1 / Gyk2 + Svyk;
    return CombFysim0;
}

double MagicFormula::PrivateCombFy_00(double alpha, double parameters[]){//checked by matlab
    const double Byk_ = cos(atan(parameters[1] * (alpha - parameters[2])));
    const double temp1 = Byk*k_s - Eyk*(Byk*k_s - atan(Byk*k_s));
    const double temp1_ = k_s*Byk_ - Eyk*k_s*Byk_ + Eyk / (1 + Byk*k_s*Byk*k_s)*k_s*Byk_;
    const double Gyk1_ = -sin(parameters[3] * atan(temp1))*parameters[3] / (1 + temp1*temp1)*temp1_;
    const double temp2 = Byk*Shyk - Eyk*(Byk*Shyk - atan(Byk*Shyk));
    const double temp2_ = Shyk*Byk_ - Eyk*Shyk*Byk_ + Eyk / (1 + Byk*Shyk*Byk*Shyk)*Shyk*Byk_;
    const double Gyk2_ = -sin(parameters[3] * atan(temp2))*parameters[3] / (1 + temp2*temp2)*temp2_;
    return Fyo*(Gyk1_*Gyk2 - Gyk1*Gyk2_) / (Gyk2*Gyk2);
}

double MagicFormula::PrivateCombFy_01(double alpha, double parameters[]){//checked by matlab
    const double temp0 = parameters[1] * (alpha - parameters[2]);
    const double Byk_ = -parameters[0] * sin(atan(temp0)) / (1 + temp0*temp0)*(alpha - parameters[2]);
    const double temp1 = Byk*k_s - Eyk*(Byk*k_s - atan(Byk*k_s));
    const double temp1_ = k_s*Byk_ - Eyk*k_s*Byk_ + Eyk / (1 + Byk*k_s*Byk*k_s)*k_s*Byk_;
    const double Gyk1_ = -sin(parameters[3] * atan(temp1))*parameters[3] / (1 + temp1*temp1)*temp1_;
    const double temp2 = Byk*Shyk - Eyk*(Byk*Shyk - atan(Byk*Shyk));
    const double temp2_ = Shyk*Byk_ - Eyk*Shyk*Byk_ + Eyk / (1 + Byk*Shyk*Byk*Shyk)*Shyk*Byk_;
    const double Gyk2_ = -sin(parameters[3] * atan(temp2))*parameters[3] / (1 + temp2*temp2)*temp2_;
    return Fyo*(Gyk1_*Gyk2 - Gyk1*Gyk2_) / (Gyk2*Gyk2);
}

double MagicFormula::PrivateCombFy_02(double alpha, double parameters[]){//checked by matlab
    const double temp0 = parameters[1] * (alpha - parameters[2]);
    const double Byk_ = parameters[0] * sin(atan(temp0)) / (1 + temp0*temp0)*parameters[1];
    const double temp1 = Byk*k_s - Eyk*(Byk*k_s - atan(Byk*k_s));
    const double temp1_ = k_s*Byk_ - Eyk*k_s*Byk_ + Eyk / (1 + Byk*k_s*Byk*k_s)*k_s*Byk_;
    const double Gyk1_ = -sin(parameters[3] * atan(temp1))*parameters[3] / (1 + temp1*temp1)*temp1_;
    const double temp2 = Byk*Shyk - Eyk*(Byk*Shyk - atan(Byk*Shyk));
    const double temp2_ = Shyk*Byk_ - Eyk*Shyk*Byk_ + Eyk / (1 + Byk*Shyk*Byk*Shyk)*Shyk*Byk_;
    const double Gyk2_ = -sin(parameters[3] * atan(temp2))*parameters[3] / (1 + temp2*temp2)*temp2_;
    return Fyo*(Gyk1_*Gyk2 - Gyk1*Gyk2_) / (Gyk2*Gyk2);
}

double MagicFormula::PrivateCombFy_03(double parameters[]){//checked by matlab
    const double temp1 = atan(Byk*k_s - Eyk*(Byk*k_s - atan(Byk*k_s)));
    const double Gyk1_ = -sin(parameters[3] * temp1)*temp1;
    const double temp2 = atan(Byk*Shyk - Eyk*(Byk*Shyk - atan(Byk*Shyk)));
    const double Gyk2_ = -sin(parameters[3] * temp2)*temp2;
    return Fyo*(Gyk1_*Gyk2 - Gyk1*Gyk2_) / (Gyk2*Gyk2);
}

double MagicFormula::PrivateCombFy_04(double parameters[]){//checked by matlab
    const double Eyk_ = 1;
    const double temp1 = Byk*k_s - Eyk*(Byk*k_s - atan(Byk*k_s));
    const double temp1_ = -(Byk*k_s - atan(Byk*k_s))*Eyk_;
    const double Gyk1_ = -sin(parameters[3] * atan(temp1))*parameters[3] / (1 + temp1*temp1)*temp1_;
    const double temp2 = Byk*Shyk - Eyk*(Byk*Shyk - atan(Byk*Shyk));
    const double temp2_ = -(Byk*Shyk - atan(Byk*Shyk))*Eyk_;
    const double Gyk2_ = -sin(parameters[3] * atan(temp2))*parameters[3] / (1 + temp2*temp2)*temp2_;
    return Fyo*(Gyk1_*Gyk2 - Gyk1*Gyk2_) / (Gyk2*Gyk2);
}

double MagicFormula::PrivateCombFy_05(double parameters[]){//checked by matlab
    const double Eyk_ = dfz;
    const double temp1 = Byk*k_s - Eyk*(Byk*k_s - atan(Byk*k_s));
    const double temp1_ = -(Byk*k_s - atan(Byk*k_s))*Eyk_;
    const double Gyk1_ = -sin(parameters[3] * atan(temp1))*parameters[3] / (1 + temp1*temp1)*temp1_;
    const double temp2 = Byk*Shyk - Eyk*(Byk*Shyk - atan(Byk*Shyk));
    const double temp2_ = -(Byk*Shyk - atan(Byk*Shyk))*Eyk_;
    const double Gyk2_ = -sin(parameters[3] * atan(temp2))*parameters[3] / (1 + temp2*temp2)*temp2_;
    return Fyo*(Gyk1_*Gyk2 - Gyk1*Gyk2_) / (Gyk2*Gyk2);
}

double MagicFormula::PrivateCombFy_06(double parameters[]){//checked by matlab
    const double Shyk_ = 1;
    const double k_s_ = 1;
    const double temp1 = Byk*k_s - Eyk*(Byk*k_s - atan(Byk*k_s));
    const double temp1_ = Byk*k_s_ - Eyk*Byk*k_s_ + Eyk / (1 + Byk*k_s*Byk*k_s)*Byk*k_s_;
    const double Gyk1_ = -sin(parameters[3] * atan(temp1))*parameters[3] / (1 + temp1*temp1)*temp1_;
    const double temp2 = Byk*Shyk - Eyk*(Byk*Shyk - atan(Byk*Shyk));
    const double temp2_ = Byk*Shyk_ - Eyk*Byk*Shyk_ + Eyk / (1 + Byk*Shyk*Byk*Shyk)*Byk*Shyk_;
    const double Gyk2_ = -sin(parameters[3] * atan(temp2))*parameters[3] / (1 + temp2*temp2)*temp2_;
    return Fyo*(Gyk1_*Gyk2 - Gyk1*Gyk2_) / (Gyk2*Gyk2);
}

double MagicFormula::PrivateCombFy_07(double parameters[]){//checked by matlab
    const double Shyk_ = dfz;
    const double k_s_ = dfz;
    const double temp1 = Byk*k_s - Eyk*(Byk*k_s - atan(Byk*k_s));
    const double temp1_ = Byk*k_s_ - Eyk*Byk*k_s_ + Eyk / (1 + Byk*k_s*Byk*k_s)*Byk*k_s_;
    const double Gyk1_ = -sin(parameters[3] * atan(temp1))*parameters[3] / (1 + temp1*temp1)*temp1_;
    const double temp2 = Byk*Shyk - Eyk*(Byk*Shyk - atan(Byk*Shyk));
    const double temp2_ = Byk*Shyk_ - Eyk*Byk*Shyk_ + Eyk / (1 + Byk*Shyk*Byk*Shyk)*Byk*Shyk_;
    const double Gyk2_ = -sin(parameters[3] * atan(temp2))*parameters[3] / (1 + temp2*temp2)*temp2_;
    return Fyo*(Gyk1_*Gyk2 - Gyk1*Gyk2_) / (Gyk2*Gyk2);
}

double MagicFormula::PrivateCombFy_08(double Fz, double k, double alpha, double parameters[]){//checked by matlab
    return uy*Fz*cos(atan(parameters[11] * alpha))*sin(parameters[12] * atan(parameters[13] * k));
}

double MagicFormula::PrivateCombFy_09(double Fz, double k, double alpha, double parameters[]){//checked by matlab
    return uy*Fz*dfz*cos(atan(parameters[11] * alpha))*sin(parameters[12] * atan(parameters[13] * k));
}

double MagicFormula::PrivateCombFy_10(double Fz, double k, double alpha, double r, double parameters[]){//checked by matlab
    return uy*Fz*r*cos(atan(parameters[11] * alpha))*sin(parameters[12] * atan(parameters[13] * k));
}

double MagicFormula::PrivateCombFy_11(double Fz, double k, double alpha, double r, double parameters[]){//checked by matlab
    const double temp_ = -sin(atan(parameters[11] * alpha)) / (1 + parameters[11] * alpha*parameters[11] * alpha)*alpha;
    return uy*Fz*(parameters[8] + parameters[9] * dfz + parameters[10] * r)*sin(parameters[12] * atan(parameters[13] * k))*temp_;
}

double MagicFormula::PrivateCombFy_12(double Fz, double k, double alpha, double r, double parameters[]){//checked by matlab
    const double temp_ = cos(parameters[12] * atan(parameters[13] * k))*atan(parameters[13] * k);
    return uy*Fz*(parameters[8] + parameters[9] * dfz + parameters[10] * r)*cos(atan(parameters[11] * alpha))*temp_;
}

double MagicFormula::PrivateCombFy_13(double Fz, double k, double alpha, double r, double parameters[]){//checked by matlab
    const double temp_ = cos(parameters[12] * atan(parameters[13] * k))*parameters[12] / (1 + parameters[13] * k*parameters[13] * k)*k;
    return uy*Fz*(parameters[8] + parameters[9] * dfz + parameters[10] * r)*cos(atan(parameters[11] * alpha))*temp_;
}


double MagicFormula::PrivateCombMz(double Fz, double k, double alpha, double r, double parameters[]){//联合工况回正力矩
    /*与matlab源代码的系数对应关系表:
    Ssz1=parameters[0];
    Ssz2=parameters[1];
    Ssz3=parameters[2];
    Ssz4=parameters[3];
    */
    /*matlab源程序方程代码:
    Fzo=4000;
    dfz=(Fz-Fzo)./Fzo;
    CFk=Fz.*(Pkx1+Pkx2*dfz).*exp(Pkx3*dfz);

    Sht=qhz1+qhz2*dfz+(qhz3+qhz4*dfz).*r;
    alpha_t=tan(alpha)+Sht;
    Bt=(qbz1+qbz2*dfz+qbz3*dfz.^2).*(1+qbz4*r+qbz5*fabs(r));
    Ct=qcz1;
    Dt=Fz.*(qdz1+qdz2*dfz).*(1+qdz3*r+qdz4*r.^2).*(Ro./Fzo);
    Et=(qez1+qez2*dfz+qez3*dfz.^2).*(1+(qez4+qez5*r).*((2/pi)*atan(Bt.*Ct.*alpha_t)));
    Shy=(Phy1+Phy2*dfz)+Phy3*r-1;
    Svy=Fz.*((Pvy1+Pvy2*dfz)+(Pvy3+Pvy4*dfz).*r);
    CFyo=Pky1*Fzo.*sin(2*atan(Fz./(Pky2*Fzo)));
    CFy=CFyo.*(1-Pky3*fabs(r));
    Shf=Shy+Svy./CFy;
    alpha_r=tan(alpha)+Shf;
    Cy=Pcy1;
    uy=(Pdy1+Pdy2*dfz).*(1-Pdy3*r.^2);
    Dy=uy.*Fz;
    By=CFy./(Cy.*Dy);
    Br=qbz9+qbz10.*By.*Cy;
    Dr=Fz.*((qdz6+qdz7*dfz)+(qdz8+qdz9*dfz).*r).*Ro-1;


    Dvyk=uy.*Fz.*(Rvy1+Rvy2.*dfz+Rvy3.*r).*cos(atan(Rvy4.*alpha));
    Svyk=Dvyk.*sin(Rvy5.*atan(Rvy6.*k));
    CombFysim0=PrivateCombFy(Fz,k,alpha,r,parametersCombinedFy,parametersFy);

    alpha_teq=atan(sqrt(tan(alpha_t).^2+(CFk./CFy).^2.*k.^2.*sign(alpha_t)));
    t=Dt.*cos(Ct.*atan(Bt.*alpha_teq-Et.*(Bt.*alpha_teq-atan(Bt.*alpha_teq)))).*cos(alpha);
    CombFysim=CombFysim0-Svyk;
    alpha_req=atan(sqrt(tan(alpha_r).^2+(CFk./CFy).^2.*k.^2.*sign(alpha_r)));
    Mzr=Dr.*cos(atan(Br.*alpha_req)).*cos(alpha);
    s=Ro.*(Ssz1+Ssz2.*(CombFysim0./Fzo)+(Ssz3+Ssz4.*dfz).*r);
    CombFxsim0=PrivateCombFx(Fz,k,alpha,r,parametersCombinedFx,parametersFx);
    Mzsim=-t.*CombFysim+Mzr+s.*CombFxsim0;
    */

    PrivateMz(Fz, alpha, r, parametersPrivateMz);//计算PrivateMz的中间数据

    double temp = tan(alpha_t);
    alpha_teq = atan(sqrt(temp*temp + (CFx / CFy)*(CFx / CFy)*k*k*sign(alpha_t)));
    tMz = Dt*cos(Ct*atan(Bt*alpha_teq - Et*(Bt*alpha_teq - atan(Bt*alpha_teq))))*cos(alpha);
    temp = tan(alpha_r);
    alpha_req = atan(sqrt(temp*temp + (CFx / CFy)*(CFx / CFy)*k*k*sign(alpha_r)));
    MzrMz = Dr*cos(atan(Br*alpha_req))*cos(alpha);
    s = Ro*(parameters[0] + parameters[1] * (CombFysim0 / Fzo) + (parameters[2] + parameters[3] * dfz)*r);
    return -tMz*(CombFysim0 - Svyk) + MzrMz + s*CombFxsim0;
}

double MagicFormula::PrivateCombMz_00(){
    const double s_ = Ro;
    return s_*CombFxsim0;
}

double MagicFormula::PrivateCombMz_01(){
    const double s_ = Ro*(CombFysim0 / Fzo);
    return s_*CombFxsim0;
}

double MagicFormula::PrivateCombMz_02(double r){
    const double s_ = Ro*r;
    return s_*CombFxsim0;
}

double MagicFormula::PrivateCombMz_03(double r){
    const double s_ = Ro*dfz*r;
    return s_*CombFxsim0;
}

#undef Fzo
#undef dfz

#undef Shx
#undef Svx
#undef kx
#undef Cx
#undef Dx
#undef Ex
#undef CFx
#undef Bx
#undef Fxo

#undef Shy
#undef Svy
#undef alpha_y
#undef Cy
#undef uy
#undef Dy
#undef Ey
#undef CFy
#undef By
#undef Fyo

#undef Ro
#undef Cr
#undef alpha_t
#undef Bt
#undef Ct
#undef Dt
#undef Et
#undef t
#undef alpha_r
#undef Br
#undef Dr
#undef Mzr

#undef alpha_s
#undef Bxa
#undef Exa
#undef Gxa
#undef Gxa2
#undef CombFxsim0

#undef Shyk
#undef k_s
#undef Byk
#undef Eyk
#undef Svyk
#undef Gyk1
#undef Gyk2
#undef CombFysim0

#undef alpha_teq
#undef tMz
#undef alpha_req
#undef MzrMz
#undef s
